

# to manipulate data
library(tidyverse)

# to estimate penalized likelihood models
library(brglm2)

# to cluster standard errors
library(lmtest)
library(sandwich)

# to generate results tables
library(stargazer)



# load data
load("full_data_glm.RData")


########### PERIOD: 1995-2019
# SHORT TERM
brglm1 <- glm(short_term ~ imf_program_lag + 
                previous_short_term + previous_long_term + 
                log_gdp_pc_constant_lag + gdp_growth_lag + resource_rents_lag + 
                working_age_pop_lag + field_discovery_lag + oil_price_lag + crisis_lag + wb_extractive_project_lag +
                polity2_lag + left_executive_lag + 
                election_lag + nationalization_oc_lag + 
                any_war_lag + 
                eiti_member_lag +
                as.factor(iso3c) + poly(t, 3), 
              family = binomial(link="logit"), data = subset(glm_df, year>1994), method = "brglmFit")

brglm1_r <- brglm1 %>% 
  coeftest(vcovHC(brglm1, type = 'HC0', cluster = c('iso3c', 'Summary')))


brglm2 <- glm(short_term ~ imf_program_lag + resource_conditionality_lag +
                previous_short_term + previous_long_term + 
                log_gdp_pc_constant_lag + gdp_growth_lag + resource_rents_lag + 
                working_age_pop_lag + field_discovery_lag + oil_price_lag + crisis_lag + wb_extractive_project_lag +
                polity2_lag + left_executive_lag + 
                election_lag + nationalization_oc_lag + 
                any_war_lag + 
                eiti_member_lag +
                as.factor(iso3c) + poly(t, 3), 
              family = binomial(link="logit"), data = subset(glm_df, year>1994), method = "brglmFit")

brglm2_r <- brglm2 %>% 
  coeftest(vcovHC(brglm2, type = 'HC0', cluster = c('iso3c', 'Summary')))


brglm3 <- glm(short_term ~ imf_program_lag + resource_conditionality_dict_tfidf_lag +
                previous_short_term + previous_long_term + 
                log_gdp_pc_constant_lag + gdp_growth_lag + resource_rents_lag +
                working_age_pop_lag + field_discovery_lag + oil_price_lag + crisis_lag + wb_extractive_project_lag +
                polity2_lag + left_executive_lag + 
                election_lag + nationalization_oc_lag + 
                any_war_lag + 
                eiti_member_lag +
                as.factor(iso3c) + poly(t, 3), 
              family = binomial(link="logit"), data = subset(glm_df, year>1994), method = "brglmFit")

brglm3_r <- brglm3 %>% 
  coeftest(vcovHC(brglm3, type = 'HC0', cluster = c('iso3c', 'Summary')))


# LONG TERM
brglm4 <- glm(long_term ~ imf_program_lag + 
                previous_short_term + previous_long_term + 
                log_gdp_pc_constant_lag + gdp_growth_lag + resource_rents_lag + 
                working_age_pop_lag + field_discovery_lag + oil_price_lag + crisis_lag + wb_extractive_project_lag +
                polity2_lag + left_executive_lag + 
                election_lag + nationalization_oc_lag +
                any_war_lag + 
                eiti_member_lag +
                as.factor(iso3c) + poly(t, 3), 
              family = binomial(link="logit"), data = subset(glm_df, year>1994), method = "brglmFit")

brglm4_r <- brglm4 %>% 
  coeftest(vcovHC(brglm4, type = 'HC0', cluster = c('iso3c', 'Summary')))


brglm5 <- glm(long_term ~ imf_program_lag + resource_conditionality_lag +
                previous_short_term + previous_long_term + 
                log_gdp_pc_constant_lag + gdp_growth_lag + resource_rents_lag + 
                working_age_pop_lag + field_discovery_lag + oil_price_lag + crisis_lag + wb_extractive_project_lag +
                polity2_lag + left_executive_lag + 
                election_lag + nationalization_oc_lag +
                any_war_lag + 
                eiti_member_lag +
                as.factor(iso3c) + poly(t, 3), 
              family = binomial(link="logit"), data = subset(glm_df, year>1994), method = "brglmFit")

brglm5_r <- brglm5 %>% 
  coeftest(vcovHC(brglm5, type = 'HC0', cluster = c('iso3c', 'Summary')))


brglm6 <- glm(long_term ~ imf_program_lag + resource_conditionality_dict_tfidf_lag +
                previous_short_term + previous_long_term + 
                log_gdp_pc_constant_lag + gdp_growth_lag + resource_rents_lag + 
                working_age_pop_lag + field_discovery_lag + oil_price_lag + crisis_lag + wb_extractive_project_lag +
                polity2_lag + left_executive_lag + 
                election_lag + nationalization_oc_lag + 
                any_war_lag + 
                eiti_member_lag +
                as.factor(iso3c) + poly(t, 3), 
              family = binomial(link="logit"), data = subset(glm_df, year>1994), method = "brglmFit")

brglm6_r <- brglm6 %>% 
  coeftest(vcovHC(brglm6, type = 'HC0', cluster = c('iso3c', 'Summary')))

# TABLE G.1 (with clustered standard errors)
stargazer(brglm1_r,brglm2_r,brglm3_r,brglm4_r,brglm5_r,brglm6_r, omit = c("iso3c"), #type = "text",
          covariate.labels = c("Program Participation = 1",
                               "Resource Conditionality (%)",
                               "Resource Conditionality (TF--IDF)",
                               "Previous Short-Term Policy",
                               "Previous Long-Term Policy",
                               "GDP per Capita (Log)",
                               "GDP Growth (%)",
                               "Resource Rents (% GDP)",
                               "Working Age Population (%)",
                               "Field Discovery = 1",
                               "Oil Price (USD)",
                               "Crisis = 1",
                               "WB Extractive Project",
                               "Democracy (Polity2)",
                               "Left Executive = 1",
                               "Election Year = 1",
                               "Oil Company Nationalization = 1",
                               "War = 1",
                               "EITI Member = 1"),
          no.space = T)

# retrieve n, ll etc
stargazer(brglm1,brglm2,brglm3,brglm4,brglm5,brglm6, omit = c("iso3c"), #type = "text",
          dep.var.labels = c("Short-Term Policy","Long-Term Policy"),
          covariate.labels = c("Program Participation = 1",
                               "Resource Conditionality (%)",
                               "Resource Conditionality (TF--IDF)",
                               "Previous Short-Term Policy",
                               "Previous Long-Term Policy",
                               "GDP per Capita (Log)",
                               "GDP Growth (%)",
                               "Resource Rents (% GDP)",
                               "Working Age Population (%)",
                               "Field Discovery = 1",
                               "Oil Price (USD)",
                               "Crisis = 1",
                               "WB Extractive Project",
                               "Democracy (Polity2)",
                               "Left Executive = 1",
                               "Election Year = 1",
                               "Oil Company Nationalization = 1",
                               "War = 1",
                               "EITI Member = 1"),
          no.space = T)


########### CONTROLLING FOR EXTREME GROWTH
glm_df <- glm_df %>%
  mutate(extreme_expansion = ifelse(gdp_growth > 10, 1, 0),
         extreme_contraction = ifelse(gdp_growth < (-10), 1, 0))

# SHORT TERM
brglm7 <- glm(short_term ~ imf_program_lag + 
                previous_short_term + previous_long_term + 
                log_gdp_pc_constant_lag + gdp_growth_lag + extreme_expansion + extreme_contraction + 
                resource_rents_lag + 
                working_age_pop_lag + field_discovery_lag + oil_price_lag + crisis_lag + wb_extractive_project_lag +
                polity2_lag + left_executive_lag + 
                election_lag + nationalization_oc_lag + 
                any_war_lag + 
                eiti_member_lag +
                as.factor(iso3c) + poly(t, 3), 
              family = binomial(link="logit"), data = glm_df, method = "brglmFit")

brglm7_r <- brglm7 %>% 
  coeftest(vcovHC(brglm7, type = 'HC0', cluster = c('iso3c', 'Summary')))


brglm8 <- glm(short_term ~ imf_program_lag + resource_conditionality_lag +
                previous_short_term + previous_long_term + 
                log_gdp_pc_constant_lag + gdp_growth_lag + extreme_expansion + extreme_contraction + 
                resource_rents_lag + 
                working_age_pop_lag + field_discovery_lag + oil_price_lag + crisis_lag + wb_extractive_project_lag +
                polity2_lag + left_executive_lag + 
                election_lag + nationalization_oc_lag + 
                any_war_lag + 
                eiti_member_lag +
                as.factor(iso3c) + poly(t, 3), 
              family = binomial(link="logit"), data = glm_df, method = "brglmFit")

brglm8_r <- brglm8 %>% 
  coeftest(vcovHC(brglm8, type = 'HC0', cluster = c('iso3c', 'Summary')))


brglm9 <- glm(short_term ~ imf_program_lag + resource_conditionality_dict_tfidf_lag +
                previous_short_term + previous_long_term + 
                log_gdp_pc_constant_lag + gdp_growth_lag + extreme_expansion + extreme_contraction + 
                resource_rents_lag +
                working_age_pop_lag + field_discovery_lag + oil_price_lag + crisis_lag + wb_extractive_project_lag +
                polity2_lag + left_executive_lag + 
                election_lag + nationalization_oc_lag + 
                any_war_lag + 
                eiti_member_lag +
                as.factor(iso3c) + poly(t, 3), 
              family = binomial(link="logit"), data = glm_df, method = "brglmFit")

brglm9_r <- brglm9 %>% 
  coeftest(vcovHC(brglm9, type = 'HC0', cluster = c('iso3c', 'Summary')))


# LONG TERM
brglm10 <- glm(long_term ~ imf_program_lag + 
                previous_short_term + previous_long_term + 
                log_gdp_pc_constant_lag + gdp_growth_lag + extreme_expansion + extreme_contraction + 
                resource_rents_lag + 
                working_age_pop_lag + field_discovery_lag + oil_price_lag + crisis_lag + wb_extractive_project_lag +
                polity2_lag + left_executive_lag + 
                election_lag + nationalization_oc_lag + 
                any_war_lag + 
                eiti_member_lag +
                as.factor(iso3c) + poly(t, 3), 
              family = binomial(link="logit"), data = glm_df, method = "brglmFit")

brglm10_r <- brglm10 %>% 
  coeftest(vcovHC(brglm10, type = 'HC0', cluster = c('iso3c', 'Summary')))


brglm11 <- glm(long_term ~ imf_program_lag + resource_conditionality_lag +
                previous_short_term + previous_long_term + 
                log_gdp_pc_constant_lag + gdp_growth_lag + extreme_expansion + extreme_contraction + 
                resource_rents_lag + 
                working_age_pop_lag + field_discovery_lag + oil_price_lag + crisis_lag + wb_extractive_project_lag +
                polity2_lag + left_executive_lag + 
                election_lag + nationalization_oc_lag + 
                any_war_lag + 
                eiti_member_lag +
                as.factor(iso3c) + poly(t, 3), 
              family = binomial(link="logit"), data = glm_df, method = "brglmFit")

brglm11_r <- brglm11 %>% 
  coeftest(vcovHC(brglm11, type = 'HC0', cluster = c('iso3c', 'Summary')))


brglm12 <- glm(long_term ~ imf_program_lag + resource_conditionality_dict_tfidf_lag +
                previous_short_term + previous_long_term + 
                log_gdp_pc_constant_lag + gdp_growth_lag + extreme_expansion + extreme_contraction + 
                resource_rents_lag + 
                working_age_pop_lag + field_discovery_lag + oil_price_lag + crisis_lag + wb_extractive_project_lag +
                polity2_lag + left_executive_lag + 
                election_lag + nationalization_oc_lag + 
                any_war_lag + 
                eiti_member_lag +
                as.factor(iso3c) + poly(t, 3), 
              family = binomial(link="logit"), data = glm_df, method = "brglmFit")

brglm12_r <- brglm12 %>% 
  coeftest(vcovHC(brglm12, type = 'HC0', cluster = c('iso3c', 'Summary')))



# TABLE G.2 (with clustered standard errors)
stargazer(brglm7_r,brglm8_r,brglm9_r,brglm10_r,brglm11_r,brglm12_r, omit = c("iso3c"), #type = "text",
          covariate.labels = c("Program Participation = 1",
                               "Resource Conditionality (%)",
                               "Resource Conditionality (TF--IDF)",
                               "Previous Short-Term Policy",
                               "Previous Long-Term Policy",
                               "GDP per Capita (Log)",
                               "GDP Growth (%)",
                               "Extreme Expansion = 1",
                               "Extreme Contraction = 1",
                               "Resource Rents (% GDP)",
                               "Working Age Population (%)",
                               "Field Discovery = 1",
                               "Oil Price (USD)",
                               "Crisis = 1",
                               "WB Extractive Project",
                               "Democracy (Polity2)",
                               "Left Executive = 1",
                               "Election Year = 1",
                               "Oil Company Nationalization = 1",
                               "War = 1",
                               "EITI Member = 1"),
          no.space = T)

# retrieve n, ll etc
stargazer(brglm7,brglm8,brglm9,brglm10,brglm11,brglm12, omit = c("iso3c"), #type = "text",
          dep.var.labels = c("Short-Term Policy","Long-Term Policy"),
          covariate.labels = c("Program Participation = 1",
                               "Resource Conditionality (%)",
                               "Resource Conditionality (TF--IDF)",
                               "Previous Short-Term Policy",
                               "Previous Long-Term Policy",
                               "GDP per Capita (Log)",
                               "GDP Growth (%)",
                               "Extreme Expansion = 1",
                               "Extreme Contraction = 1",
                               "Resource Rents (% GDP)",
                               "Working Age Population (%)",
                               "Field Discovery = 1",
                               "Oil Price (USD)",
                               "Crisis = 1",
                               "WB Extractive Project",
                               "Democracy (Polity2)",
                               "Left Executive = 1",
                               "Election Year = 1",
                               "Oil Company Nationalization = 1",
                               "War = 1",
                               "EITI Member = 1"),
          no.space = T)



########### CONTROLLING FOR EXTERNAL DEBT

# SHORT TERM
brglm13 <- glm(short_term ~ imf_program_lag + 
                previous_short_term + previous_long_term + 
                log_gdp_pc_constant_lag + gdp_growth_lag +
                resource_rents_lag + 
                log_external_debt_stocks_lag +
                working_age_pop_lag + field_discovery_lag + oil_price_lag + crisis_lag + wb_extractive_project_lag +
                polity2_lag + left_executive_lag + 
                election_lag + nationalization_oc_lag + 
                any_war_lag + 
                eiti_member_lag +
                as.factor(iso3c) + poly(t, 3), 
              family = binomial(link="logit"), data = glm_df, method = "brglmFit")

brglm13_r <- brglm13 %>% 
  coeftest(vcovHC(brglm13, type = 'HC0', cluster = c('iso3c', 'Summary')))


brglm14 <- glm(short_term ~ imf_program_lag + resource_conditionality_lag +
                previous_short_term + previous_long_term + 
                log_gdp_pc_constant_lag + gdp_growth_lag +
                resource_rents_lag + log_external_debt_stocks_lag +
                working_age_pop_lag + field_discovery_lag + oil_price_lag + crisis_lag + wb_extractive_project_lag +
                polity2_lag + left_executive_lag + 
                election_lag + nationalization_oc_lag + 
                any_war_lag + 
                eiti_member_lag +
                as.factor(iso3c) + poly(t, 3), 
              family = binomial(link="logit"), data = glm_df, method = "brglmFit")

brglm14_r <- brglm14 %>% 
  coeftest(vcovHC(brglm14, type = 'HC0', cluster = c('iso3c', 'Summary')))


brglm15 <- glm(short_term ~ imf_program_lag + resource_conditionality_dict_tfidf_lag +
                previous_short_term + previous_long_term + 
                log_gdp_pc_constant_lag + gdp_growth_lag +
                resource_rents_lag + log_external_debt_stocks_lag +
                working_age_pop_lag + field_discovery_lag + oil_price_lag + crisis_lag + wb_extractive_project_lag +
                polity2_lag + left_executive_lag + 
                election_lag + nationalization_oc_lag + 
                any_war_lag + 
                eiti_member_lag +
                as.factor(iso3c) + poly(t, 3), 
              family = binomial(link="logit"), data = glm_df, method = "brglmFit")

brglm15_r <- brglm15 %>% 
  coeftest(vcovHC(brglm15, type = 'HC0', cluster = c('iso3c', 'Summary')))


# LONG TERM
brglm16 <- glm(long_term ~ imf_program_lag +
                 previous_short_term + previous_long_term + 
                 log_gdp_pc_constant_lag + gdp_growth_lag +
                 resource_rents_lag + log_external_debt_stocks_lag +
                 working_age_pop_lag + field_discovery_lag + oil_price_lag + crisis_lag + wb_extractive_project_lag +
                 polity2_lag + left_executive_lag + 
                 election_lag + nationalization_oc_lag + 
                 any_war_lag + 
                 eiti_member_lag +
                 as.factor(iso3c) + poly(t, 3), 
               family = binomial(link="logit"), data = glm_df, method = "brglmFit")

brglm16_r <- brglm16 %>% 
  coeftest(vcovHC(brglm16, type = 'HC0', cluster = c('iso3c', 'Summary')))


brglm17 <- glm(long_term ~ imf_program_lag + resource_conditionality_lag +
                 previous_short_term + previous_long_term + 
                 log_gdp_pc_constant_lag + gdp_growth_lag +
                 resource_rents_lag + log_external_debt_stocks_lag +
                 working_age_pop_lag + field_discovery_lag + oil_price_lag + crisis_lag + wb_extractive_project_lag +
                 polity2_lag + left_executive_lag + 
                 election_lag + nationalization_oc_lag + 
                 any_war_lag + 
                 eiti_member_lag +
                 as.factor(iso3c) + poly(t, 3), 
               family = binomial(link="logit"), data = glm_df, method = "brglmFit")

brglm17_r <- brglm17 %>% 
  coeftest(vcovHC(brglm17, type = 'HC0', cluster = c('iso3c', 'Summary')))


brglm18 <- glm(long_term ~ imf_program_lag + resource_conditionality_dict_tfidf_lag +
                 previous_short_term + previous_long_term + 
                 log_gdp_pc_constant_lag + gdp_growth_lag +
                 resource_rents_lag + log_external_debt_stocks_lag +
                 working_age_pop_lag + field_discovery_lag + oil_price_lag + crisis_lag + wb_extractive_project_lag +
                 polity2_lag + left_executive_lag + 
                 election_lag + nationalization_oc_lag + 
                 any_war_lag + 
                 eiti_member_lag +
                 as.factor(iso3c) + poly(t, 3), 
               family = binomial(link="logit"), data = glm_df, method = "brglmFit")

brglm18_r <- brglm18 %>% 
  coeftest(vcovHC(brglm18, type = 'HC0', cluster = c('iso3c', 'Summary')))

# TABLE G.3 (with clustered standard errors)
stargazer(brglm13_r,brglm14_r,brglm15_r,brglm16_r,brglm17_r,brglm18_r, omit = c("iso3c"), #type = "text",
          covariate.labels = c("Program Participation = 1",
                               "Resource Conditionality (%)",
                               "Resource Conditionality (TF--IDF)",
                               "Previous Short-Term Policy",
                               "Previous Long-Term Policy",
                               "GDP per Capita (Log)",
                               "GDP Growth (%)",
                               "Resource Rents (% GDP)",
                               "External Debt Stocks",
                               "Working Age Population (%)",
                               "Field Discovery = 1",
                               "Oil Price (USD)",
                               "Crisis = 1",
                               "WB Extractive Project",
                               "Democracy (Polity2)",
                               "Left Executive = 1",
                               "Election Year = 1",
                               "Oil Company Nationalization = 1",
                               "War = 1",
                               "EITI Member = 1"),
          no.space = T)



# retrieve n, ll etc
stargazer(brglm13,brglm14,brglm15,brglm16,brglm17,brglm18, omit = c("iso3c"), #type = "text",
          dep.var.labels = c("Short-Term Policy","Long-Term Policy"),
          covariate.labels = c("Program Participation = 1",
                               "Resource Conditionality (%)",
                               "Resource Conditionality (TF--IDF)",
                               "Previous Short-Term Policy",
                               "Previous Long-Term Policy",
                               "GDP per Capita (Log)",
                               "GDP Growth (%)",
                               "Resource Rents (% GDP)",
                               "External Debt Stocks",
                               "Working Age Population (%)",
                               "Field Discovery = 1",
                               "Oil Price (USD)",
                               "Crisis = 1",
                               "WB Extractive Project",
                               "Democracy (Polity2)",
                               "Left Executive = 1",
                               "Election Year = 1",
                               "Oil Company Nationalization = 1",
                               "War = 1",
                               "EITI Member = 1"),
          no.space = T)
